home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
FishMarket 1.0
/
FishMarket v1.0.iso
/
fishies
/
201-225
/
disk_201
/
draco
/
src
/
complex.d
next >
Wrap
Text File
|
1992-05-06
|
9KB
|
405 lines
#drinc:float.g
#drinc:util.g
type Complex_t = struct {
float cmpl_real;
float cmpl_imag;
};
uint STACK_HEIGHT = 20;
[STACK_HEIGHT] Complex_t Stack;
*Complex_t StackTop;
/*
* cmplInitialize - set up the world.
*/
proc cmplInitialize()void:
StackTop := &Stack[STACK_HEIGHT];
corp;
/*
* _cmplpsh - calls generated by the compiler to push complex values onto our
* internal stack. The address of the complex value will be on the normal
* stack as a parameter.
*/
proc _cmplpsh(*Complex_t n)void:
StackTop := StackTop - sizeof(Complex_t);
StackTop* := n*;
corp;
/*
* _cmplpop - calls generated by the compiler to pop complex values from our
* internal stack into a destination variable. The address of the
* destination variable will be on the normal stack as a parameter.
*/
proc _cmplpop(*Complex_t n)void:
n* := StackTop*;
StackTop := StackTop + sizeof(Complex_t);
corp;
/*
* _cmplput - calls generated by the compiler to write a complex value in
* text mode.
*/
proc _cmplput()void:
write(*; '(', StackTop*.cmpl_real, ',', StackTop*.cmpl_imag, ')');
StackTop := StackTop + sizeof(Complex_t);
corp;
/*
* _cmplget - calls generated by the compiler to read a complex value in
* text mode. Various error conditions can occur. An option here would
* be to make the parentheses and the comma required. The value read in
* is left on our internal stack.
*/
proc _cmplget()void:
extern
_d_channelSkip()void,
_d_channelGetChar()char,
_d_channelUnGetChar(char ch)void,
_d_channelHadError()bool,
_d_channelError(ushort err)void;
register *Complex_t p;
register char ch;
bool gotParen;
p := StackTop - sizeof(Complex_t);
StackTop := p;
_d_channelSkip();
ch := _d_channelGetChar();
if ch = '(' then
gotParen := true;
else
_d_channelUnGetChar(ch);
gotParen := false;
fi;
read(*; p*.cmpl_real);
if not _d_channelHadError() then
_d_channelSkip();
ch := _d_channelGetChar();
if ch ~= ',' then
_d_channelUnGetChar(ch);
fi;
read(*; p*.cmpl_imag);
if not _d_channelHadError() then
if gotParen then
_d_channelSkip();
ch := _d_channelGetChar();
if ch ~= ')' then
_d_channelUnGetChar(ch);
_d_channelError(CH_BADCHAR);
fi;
fi;
fi;
fi;
corp;
/*
* _cmplneg - compiler generated for the '-' unary operator.
*/
proc _cmplneg()void:
register *Complex_t p;
p := StackTop;
p*.cmpl_real := - p*.cmpl_real;
p*.cmpl_imag := - p*.cmpl_imag;
corp;
/*
* cmpladd - compiler generated for the '+' binary operator.
*/
proc _cmpladd()void:
register *Complex_t stackTop;
stackTop := StackTop;
(stackTop + sizeof(Complex_t))*.cmpl_real :=
(stackTop + sizeof(Complex_t))*.cmpl_real + stackTop*.cmpl_real;
(stackTop + sizeof(Complex_t))*.cmpl_imag :=
(stackTop + sizeof(Complex_t))*.cmpl_imag + stackTop*.cmpl_imag;
StackTop := stackTop + sizeof(Complex_t);
corp;
/*
* _cmplsub - compler generated for the '-' binary operator.
*/
proc _cmplsub()void:
register *Complex_t stackTop;
stackTop := StackTop;
(stackTop + sizeof(Complex_t))*.cmpl_real :=
(stackTop + sizeof(Complex_t))*.cmpl_real - stackTop*.cmpl_real;
(stackTop + sizeof(Complex_t))*.cmpl_imag :=
(stackTop + sizeof(Complex_t))*.cmpl_imag - stackTop*.cmpl_imag;
StackTop := stackTop + sizeof(Complex_t);
corp;
/*
* _cmplmul - compiler generated for the '*' binary operator.
*/
proc _cmplmul()void:
register *Complex_t l, r;
float temp;
r := StackTop;
l := r + sizeof(Complex_t);
temp := l*.cmpl_real;
l*.cmpl_real := temp * r*.cmpl_real - l*.cmpl_imag * r*.cmpl_imag;
l*.cmpl_imag := l*.cmpl_imag * r*.cmpl_real + temp * r*.cmpl_imag;
StackTop := l;
corp;
/*
* _cmpldiv - compiler generated for the '/' binary operator.
*/
proc _cmpldiv()void:
register *Complex_t l, r;
float t1, t2;
r := StackTop;
l := r + sizeof(Complex_t);
t1 := l*.cmpl_real;
t2 := 1.0 / (r*.cmpl_real * r*.cmpl_real + r*.cmpl_imag * r*.cmpl_imag);
l*.cmpl_real := (t1 * r*.cmpl_real + l*.cmpl_imag * r*.cmpl_imag) * t2;
l*.cmpl_imag := (l*.cmpl_imag * r*.cmpl_real - t1 * r*.cmpl_imag) * t2;
StackTop := l;
corp;
/*****************************************************************************\
* *
* The routines from here on down could have been in another source file, in *
* which case they would have had one or more 'complex' arguments or result, *
* and would not explicitly reference 'StackTop'. This would have made them *
* slower, however, since they would have been calling '_cmplpsh', '_cmplpop', *
* etc., so they are included here to allow them direct access. All of these *
* routines are called directly by the user, and not by compiler-generated *
* code. *
* *
\*****************************************************************************/
/*
* cmpl - return a complex build from two float values.
* This differs from the compiler's ability to make a complex constant
* from two float constants in that the expressions need not be
* constants, and the construction is done at run-time.
*/
proc cmpl(float re, im)void:
StackTop := StackTop - sizeof(Complex_t);
StackTop*.cmpl_real := re;
StackTop*.cmpl_imag := im;
corp;
/*
* cmplRe - return the real part of a complex number.
*/
proc cmplRe()float:
StackTop := StackTop + sizeof(Complex_t);
(StackTop - sizeof(Complex_t))*.cmpl_real
corp;
/*
* cmplIm - return the imaginary part of a complex number.
*/
proc cmplIm()float:
StackTop := StackTop + sizeof(Complex_t);
(StackTop - sizeof(Complex_t))*.cmpl_imag
corp;
/*
* cmplConj - return the complex conjugate of a value.
*/
proc cmplConj()void:
StackTop*.cmpl_imag := - StackTop*.cmpl_imag;
corp;
/*
* cmplNorm - return the complex norm (a float value) of a complex value.
*/
proc cmplNorm()float:
register *Complex_t p;
p := StackTop;
StackTop := p + sizeof(Complex_t);
p*.cmpl_real * p*.cmpl_real + p*.cmpl_imag * p*.cmpl_imag
corp;
/*
* cmplAbs - complex absolute value (a float value) of a complex value.
*/
proc cmplAbs()float:
sqrt(cmplNorm())
corp;
/*
* cmplSqrt - complex square root of a complex value.
*/
proc cmplSqrt()void:
register *Complex_t p;
float temp;
p := StackTop;
temp := cmplAbs();
StackTop := p;
p*.cmpl_imag := sqrt((temp - p*.cmpl_real) * 0.5);
p*.cmpl_real := sqrt((temp + p*.cmpl_real) * 0.5);
corp;
/*
* cmplArg
*/
proc cmplArg()float:
register *Complex_t p;
float temp;
p := StackTop;
StackTop := p + sizeof(Complex_t);
if p*.cmpl_real = 0.0 then
if p*.cmpl_imag > 0.0 then
HALF_PI
elif p*.cmpl_imag < 0.0 then
- HALF_PI
else
0.0
fi
else
temp := atan(p*.cmpl_imag / p*.cmpl_real);
if p*.cmpl_real > 0.0 then
temp
else
temp + PI
fi
fi
corp;
/*
* cmplExp
*/
proc cmplExp()void:
register *Complex_t p;
float temp;
p := StackTop;
temp := exp(p*.cmpl_real);
p*.cmpl_real := temp * cos(p*.cmpl_imag);
p*.cmpl_imag := temp * sin(p*.cmpl_imag);
corp;
/*
* cmplLog
*/
proc cmplLog()void:
register *Complex_t p;
float temp;
p := StackTop;
temp := log(cmplNorm()) * 0.5;
StackTop := p;
p*.cmpl_imag := cmplArg();
StackTop := p;
p*.cmpl_real := temp;
corp;
/*
* cmplSin
*/
proc cmplSin()void:
register *Complex_t p;
float temp;
p := StackTop;
temp := p*.cmpl_real;
p*.cmpl_real := sin(temp) * cosh(p*.cmpl_imag);
p*.cmpl_imag := cos(temp) * sinh(p*.cmpl_imag);
corp;
/*
* cmplCos
*/
proc cmplCos()void:
register *Complex_t p;
float temp;
p := StackTop;
temp := p*.cmpl_real;
p*.cmpl_real := cos(temp) * cosh(p*.cmpl_imag);
p*.cmpl_imag := - (sin(temp) * sinh(p*.cmpl_imag));
corp;
/*
proc cmplTan()void:corp;
proc cmplAsin()void:corp;
proc cmplAcos()void:corp;
*/
/*
* atanh - utility routine needed below.
*/
proc atanh(float n)float:
0.5 * log((1.0 + n) / (1.0 - n))
corp;
/*
* cmplAtan
*/
proc cmplAtan()void:
register *Complex_t p;
float anorm, snorm;
p := StackTop;
if p*.cmpl_imag = 0.0 then
p*.cmpl_real := atan(p*.cmpl_real);
elif p*.cmpl_real = 0.0 then
if p*.cmpl_imag <= -1.0 then
p*.cmpl_imag := atanh(1.0 / p*.cmpl_imag);
p*.cmpl_real := - HALF_PI;
elif p*.cmpl_imag >= 1.0 then
p*.cmpl_imag := atanh(1.0 / p*.cmpl_imag);
p*.cmpl_real := HALF_PI;
else
p*.cmpl_imag := atanh(p*.cmpl_imag);
fi;
else
snorm := cmplNorm();
StackTop := p;
anorm := snorm - 1.0;
anorm := sqrt(anorm * anorm + 4.0 * p*.cmpl_real * p*.cmpl_real);
p*.cmpl_real := atan(0.5 * (snorm + anorm - 1.0) / p*.cmpl_real);
p*.cmpl_imag := atanh(0.5 * (snorm - anorm + 1.0) / p*.cmpl_imag);
fi;
corp;